; this file reads in the WMAP5 OWLS constants and units

; some WMAP5 OWLS file
owls_file = '/nfs/cluster/galaxy1/rcrain/OWLS/THERMAL_FEEDBACK/L100N512/Data/Snapshots/snapshot_000/snap_000.0.hdf5'

; units
;print,'  -units'
grp = 'Units'
unitlength_in_cm         = h5ra(owls_file,grp,'UnitLength_in_cm')
unitmass_in_g            = h5ra(owls_file,grp,'UnitMass_in_g')
unitvelocity_in_cm_per_s = h5ra(owls_file,grp,'UnitVelocity_in_cm_per_s')
unitdensity_in_cgs       = h5ra(owls_file,grp,'UnitDensity_in_cgs')
unitenergy_in_cgs        = h5ra(owls_file,grp,'UnitEnergy_in_cgs')
unitpressure_in_cgs      = h5ra(owls_file,grp,'UnitPressure_in_cgs')
unittime_in_s            = h5ra(owls_file,grp,'UnitTime_in_s')
units = create_struct( $
        'unitlength_in_cm'         ,unitlength_in_cm,$
        'unitmass_in_g'            ,unitmass_in_g,$
        'unitvelocity_in_cm_per_s' ,unitvelocity_in_cm_per_s,$
        'unitdensity_in_cgs'       ,unitdensity_in_cgs,$
        'unitenergy_in_cgs'        ,unitenergy_in_cgs,$
        'unitpressure_in_cgs'      ,unitpressure_in_cgs,$
        'unittime_in_s'            ,unittime_in_s)

grp = 'Header'
omega0           = h5ra(owls_file,grp,'Omega0')
omegabaryon      = h5ra(owls_file,grp,'OmegaBaryon')
omegalambda      = h5ra(owls_file,grp,'OmegaLambda')
hubbleparam      = h5ra(owls_file,grp,'HubbleParam')

; constants
grp = 'Constants'
pi               = h5ra(owls_file,grp,'PI')
gamma            = h5ra(owls_file,grp,'GAMMA')
gravity          = h5ra(owls_file,grp,'GRAVITY')
solar_mass       = h5ra(owls_file,grp,'SOLAR_MASS')
solar_lum        = h5ra(owls_file,grp,'SOLAR_LUM')
rad_const        = h5ra(owls_file,grp,'RAD_CONST')
avogadro         = h5ra(owls_file,grp,'AVOGADRO')
boltzmann        = h5ra(owls_file,grp,'BOLTZMANN')
gas_const        = h5ra(owls_file,grp,'GAS_CONST')
c                = h5ra(owls_file,grp,'C')
planck           = h5ra(owls_file,grp,'PLANCK')
cm_per_mpc       = h5ra(owls_file,grp,'CM_PER_MPC')
protonmass       = h5ra(owls_file,grp,'PROTONMASS')
electronmass     = h5ra(owls_file,grp,'ELECTRONMASS')
electroncharge   = h5ra(owls_file,grp,'ELECTRONCHARGE')
hubble           = h5ra(owls_file,grp,'HUBBLE')
t_cmb0           = h5ra(owls_file,grp,'T_CMB0')
sec_per_megayear = h5ra(owls_file,grp,'SEC_PER_MEGAYEAR')
sec_per_year     = h5ra(owls_file,grp,'SEC_PER_YEAR')
stefan           = h5ra(owls_file,grp,'STEFAN')
thompson         = h5ra(owls_file,grp,'THOMPSON')
ev_to_erg        = h5ra(owls_file,grp,'EV_TO_ERG')
z_solar          = h5ra(owls_file,grp,'Z_Solar')

; derive useful extras
rho_crit_cgs     = 3.0d0*hubble*hubble
rho_crit_cgs    /= 8.0d0*pi*gravity
rho_b_cgs        = rho_crit_cgs * omegabaryon
rho_m_cgs        = rho_crit_cgs * omega0
rho_crit         = rho_crit_cgs / units.unitdensity_in_cgs
rho_b            = rho_crit     * omegabaryon
rho_m            = rho_crit     * omega0


constants        = create_struct( $
                   'pi'               ,pi,$
                   'gamma'            ,gamma,$
                   'gravity'          ,gravity,$
                   'solar_mass'       ,solar_mass,$
                   'solar_lum'        ,solar_lum,$
                   'rad_const'        ,rad_const,$
                   'avogadro'         ,avogadro,$
                   'boltzmann'        ,boltzmann,$
                   'gas_const'        ,gas_const,$
                   'c'                ,c,$
                   'planck'           ,planck,$
                   'cm_per_mpc'       ,cm_per_mpc,$
                   'protonmass'       ,protonmass,$
                   'electronmass'     ,electronmass,$
                   'electroncharge'   ,electroncharge,$
                   'hubble'           ,hubble,$
                   't_cmb0'           ,t_cmb0,$
                   'sec_per_megayear' ,sec_per_megayear,$
                   'sec_per_year'     ,sec_per_year,$
                   'stefan'           ,stefan,$
                   'thompson'         ,thompson,$
                   'ev_to_erg'        ,ev_to_erg,$
                   'z_solar'          ,z_solar,$
                   'rho_crit'         ,rho_crit,$
                   'rho_m'            ,rho_m,$
                   'rho_b'            ,rho_b,$
                   'rho_m_cgs'        ,rho_m_cgs,$
                   'rho_b_cgs'        ,rho_b_cgs,$
                   'rho_crit_cgs'     ,rho_crit_cgs,$
                   'omega0'           ,omega0,$
                   'omegabaryon'      ,omegabaryon,$
                   'omegalambda'      ,omegalambda,$
                   'hubbleparam'      ,hubbleparam)

;print,' -done'
